Decay Rates¶

In [1]:
%matplotlib inline
In [2]:
import os,sys
sys.path.append('/home/dpirvu/python_stuff/')
sys.path.append('./bubbles_codes/')

from plotting import *
from bubble_tools import *
from experiment import *
from celluloid import Camera
In [3]:
bubbleList, velocitesList, instantonList, tmpList, fldcritList, tcritList, encritList = [], [], [], [], [], [], []
for tmp in [0,1,2,3]:
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])

    bubbleList.append(np.load(average_file(*exp_params)))
    velocitesList.append(np.load(velocities_file(*exp_params)))
    instantonList.append(np.load(supercrit_instanton_file(*exp_params)))
    tmpList.append(tmp)
    fldcritList.append(np.load(critfield_file(*exp_params)))
    tcritList.append(np.load(crittimes_file(*exp_params)))
    encritList.append(np.load(critenerg_file(*exp_params)))
In [4]:
# standardize plots
lsl = lambda tmp: ('-' if tmp%4==0 else '-.' if tmp%4==1 else '--' if tmp%4==3 else ':')
In [5]:
free_eigenbasis_mom  = lambda la, ph0: np.asarray([norm(ph0)*(w2(la)[k]**0.25) if kk!=0. else 0. for k,kk in enumerate(klist)])
thermal_eigenbasis_mom  = lambda la, ph0, te: free_eigenbasis_mom(la, ph0) * np.sqrt(2./(np.exp(w2(la)**0.5/te)-1.))
fluct_stdev_mom = lambda la, ph0, te: np.sqrt(np.sum(np.abs(thermal_eigenbasis_mom(la, ph0, te))**2.))
In [6]:
%run './bubbles_codes/plotting.py'
%run './bubbles_codes/bubble_tools.py'
In [7]:
if True:
        tmp = 0
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])
        sigmamom = fluct_stdev_mom(lamb, phi0, temp)

        print(sigmafld, sigmamom, right_Vmax.x, m2(lamb)**0.5, temp/m2(lamb)**0.5)
        print((right_Vmax.x-np.pi)/sigmafld)

        philist = np.linspace(-0.25*np.pi, 2.25*np.pi, 10000)

        fig, ax = plt.subplots(1, 1, figsize=(3.5, 2.5))
        ax.plot(philist, V(philist, lamb)/V(np.pi, lamb), label=r'$\lambda=$'+str(lamb), color='k')

        ax.set_xlabel(r'$\varphi/\varphi_0$', labelpad = 5, fontsize=15)
        ax.set_ylabel(r'$V(\varphi)/V(\varphi_{\rm fv})$', labelpad = 5, fontsize=15)
        ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi))
        ax.xaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))

        a = ax.get_yticks().tolist()[1:-1:2]
        ax.set_yticks(a)
        a[0] = r'${:.0f}$'.format(a[0])
        a[1] = r'${:.1f}$'.format(a[1])
        a[2] = r'${:.0f}$'.format(a[2])
        ax.set_yticklabels(a)
       
        ax.grid(ls=':', color='darkgray', alpha=0.5)
      # plt.legend()
        plt.savefig('./plots/potential.pdf', dpi=500)
        plt.tight_layout()
        plt.show()

0.27428630597271203 0.03977359171321589 4.251836544352683 0.1 0.8999999999999999
4.047755453286628
In [8]:
[fluct_stdev(lambList[0], phi0List[0], tempList[ii]) for ii in [0,1,2,3]]
Out[8]:
[0.27428630597271203,
 0.30326206559910246,
 0.33081957224035113,
 0.3571216688274275]
In [9]:
[fluct_stdev_mom(lambList[0], phi0List[0], tempList[ii]) for ii in [0,1,2,3]]
Out[9]:
[0.03977359171321589,
 0.04522614119404094,
 0.050656217225139674,
 0.05606566994218313]
In [10]:
[phi0List[0]/fluct_stdev(lambList[0], phi0List[0], tempList[ii]) for ii in [0,1,2,3]]
Out[10]:
[5.090532670392863, 4.6041478970906144, 4.220619088948682, 3.909769480468524]
In [11]:
[phi0List[0]/fluct_stdev_mom(lambList[0], phi0List[0], tempList[ii]) for ii in [0,1,2,3]]
Out[11]:
[35.1052882440968, 30.872928017556298, 27.563514965790336, 24.904070584286945]
In [12]:
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
In [13]:
Vbare = lambda phi, lamb: 4.*nu* (- np.cos(phi) + 0.5 * lamb**2. * np.sin(phi)**2.)
Vdoub = lambda phi, lamb: 4.*nu* (np.cos(phi) + lamb**2. * np.cos(2.*phi) )
Vtrip = lambda phi, lamb: 4.*nu* (- np.sin(phi) - 4.*lamb**2. * np.cos(phi) * np.sin(phi) )
In [14]:
philist = np.linspace(0.5*np.pi, 1.5*np.pi, 10000)
yo = Vbare(philist, lambList[0])
xo = philist

plt.plot(philist, Vbare(philist, lambList[0]))
plt.plot(philist, Vdoub(philist, lambList[0]))
plt.plot(philist, Vtrip(philist, lambList[0]))

print(Vdoub(np.pi, lamb))

dyo = np.diff(yo,1)
dxo = np.diff(xo,1)

yfirst = dxo/dyo
xfirst = 0.5*(xo[:-1]+xo[1:])

dyfirst = np.diff(yfirst,1)
dxfirst = np.diff(xfirst,1)
ysecond = dyfirst/dxfirst

xsecond=0.5*(xfirst[:-1]+xfirst[1:])

#plt.plot(x[1:-1], xsecond)

0.01

/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:14: RuntimeWarning: divide by zero encountered in true_divide
In [15]:
if True:
    for tmp in [0,1,2,3]:
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])
        sigmamom = fluct_stdev_mom(lamb, phi0, temp)

        print(sigmafld, sigmamom, right_Vmax.x, m2(lamb)**0.5, temp/m2(lamb)**0.5)

        philist = np.linspace(0.75*np.pi, 2.25*np.pi, 10000)

        fig, ax = plt.subplots(1, 1, figsize=(5, 2.5))
       # ax.plot(philist, phi0**2.*V(philist, lamb), label=r'$\lambda=$'+str(lamb), color='k')
        ax.plot(philist, V(philist, lamb)/V(np.pi, lamb), label=r'$\lambda=$'+str(lamb), color='k')
       # ax.axvline(right_Vmax.x, ls='-', color='k', linewidth=0.5)
        ax.axvline((4.*sigmafld+np.pi), ls='-', color='k', linewidth=0.5)
        for nn in np.linspace(0, 5, 6):
            ax.plot(np.pi + nn*sigmafld, V(np.pi + nn*sigmafld, lamb)/V(np.pi, lamb), \
                    marker='o', ms=3, color='r', label=int(nn))
        ax.plot(np.pi + phi0/sigmafld, V(np.pi + phi0/sigmafld, lamb)/V(np.pi, lamb), \
                marker='*', ms=3, color='g', label=int(nn))

        ax.set_xlabel(r'$\varphi/\phi_0$', labelpad = 10)
        ax.set_ylabel(r'$V(\varphi)/V(\varphi_{\rm fv})$', labelpad = 10)
        ax.xaxis.set_major_locator(plt.MultipleLocator(np.pi))
        ax.xaxis.set_major_formatter(plt.FuncFormatter(multiple_formatter()))

        a = ax.get_xticks().tolist()[1:]
        a[2] = a[1]
        a[1] = (4.*sigmafld+np.pi)#right_Vmax.x
        ax.set_xticks(a)
        a[0] = r'$\pi$'
        a[1] = r'${:.0f}$'.format((right_Vmax.x-np.pi)/sigmafld)+r'$\, \sigma_\varphi$'
        a[2] = r'$2\pi$'
        ax.set_xticklabels(a)

        ax.grid(ls=':', color='darkgray', alpha=0.5)
        plt.legend(bbox_to_anchor=(1., 1))
       # plt.savefig('./plots/potential.pdf', dpi=500)
        plt.tight_layout()
        plt.show()

0.27428630597271203 0.03977359171321589 4.251836544352683 0.1 0.8999999999999999

0.30326206559910246 0.04522614119404094 4.251836544352683 0.1 1.0

0.33081957224035113 0.050656217225139674 4.251836544352683 0.1 1.0999999999999999

0.3571216688274275 0.05606566994218313 4.251836544352683 0.1 1.2
In [16]:
if False:
    for tmp in range(len(tempList)):
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        def fun(y, x):
            Φ, Π = y
            #dydx = [Π, -Π/x+(np.sin(Φ) + 0.5*lamb**2.*np.sin(2.*Φ))*4.*nu]
            # no friction term in thermal case
            dydx = [Π, +(np.sin(Φ) + 0.5*lamb**2.*np.sin(2.*Φ))*4.*nu]
            return dydx

        #y0 = [5.69219082899233042025, 0.] # with friction
        y0 = [4.823729806624866, 0.] # supercrit_instanton_file
        
        # as far as fortran goes:
             #4.8238 is super-critical
             #4.823729806620144 is sub-critical
            
        #y0 = [4.823729806625, 0.] # instanton_file

        x     = np.linspace(0, lenLat, 1000)
        xplot = np.linspace(np.pi*0.8, 2.*np.pi, 1000)
        sol   = odeint(fun, y0, x)

        plt.plot(x*np.sqrt(m2(lamb)), sol[:, 0], ls='-', label=r'$\phi(x)$')
        plt.axhline(np.pi, ls=':', color='darkgray')
        plt.legend(); plt.show()

        plt.plot(xplot, Vinv(xplot, lamb), label='V', linewidth='1', alpha=0.3)
        plt.plot(sol[:, 0], Vinv(sol[:, 0], lamb))
        plt.legend(); plt.show()

        ind  = np.argmin(np.abs(sol[:,0] - np.pi)); print(ind, sol[ind,0])
        inst = np.ones(nLat//2+1) * phieq
        inst[:ind+1] = sol[:ind+1, 0]
        inst = np.concatenate((inst[::-1],inst[1:-1]))
        print(np.argmax(inst))
        
       # np.save(supercrit_instanton_file(*exp_params), inst)

        plt.plot(np.arange(nLat)-nLat//2, inst); plt.show()
        print('Instanton saved.')
In [17]:
if True:
    simLists = []
    for tmp in [0,1,2,3]:
        simLists.append([])
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])
        for sim in range(minSim, maxSim):
            path2CLEANsim = clean_sim_location(*exp_params, sim)
            path2RESTsim = rest_sim_location(*exp_params, sim)
            if os.path.exists(path2RESTsim) and os.path.exists(path2CLEANsim):
                simLists[-1].append(sim)
    tp = 0# 0 for average, 1 for error
    cp = 0

    titl = [r'$\bar{\varphi}$',  r'$\bar{\Pi}$', r'$\partial_x\bar{\varphi}$']

    crit_rad = 40
    crit_thresh = right_Vmax.x + 3.*sigmafld
    win = 180

    for tmp in [0,1,2,3]:
        if tmp!=0: break

        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        dats = random.sample(simLists[tmp], 6)
        fig, ax = plt.subplots(2,3, figsize = (5.5,3.5), sharey=True, sharex=True, constrained_layout=True)
        for si, sim in enumerate(dats):
            ci,ti = np.divmod(si,3)

            path2CLEANsim = clean_sim_location(*exp_params, sim)
            fullreal, sim, tdecay, outcome = np.load(path2CLEANsim)

            bubble2measure = fullreal[0]
            nT, nN = np.shape(bubble2measure)
            tcen, xcen = find_nucleation_center(bubble2measure, phieq, crit_thresh, crit_rad)
            tcen -= 100
            tl,tr = max(0, tcen-win), min(nT-1, tcen+win)
            xl,xr = max(0, xcen-win), min(nN-1, xcen+win)
            ext = np.array([xl-xcen,xr-xcen,tl-tcen-100,tr-tcen-100])*np.sqrt(m2(lamb))*dx

            bubble2plot = bubble2measure[tl:tr,xl:xr]
            im = ax[ci,ti].imshow(bubble2plot, interpolation=None, extent=ext, origin='lower', cmap='RdBu')

            if False:
                lavs = 6
                nT, nN = np.shape(bubble2plot)
                tt = np.linspace(tl-tcen, tr-tcen, nT)*np.sqrt(m2(lamb))*dx
                xx = np.linspace(xl-xcen, xr-xcen, nN)*np.sqrt(m2(lamb))*dx
                ttt1, xxx1 = np.meshgrid(tt, xx)
                ax[ci,ti].contour(xxx1, ttt1, bubble2plot.T, levels=lavs, aspect='auto', interpolation=None, extent=ext, origin='lower', colors='k',linewidths=0.5)

            if ci==1: ax[ci,ti].set(xlabel=r'$m \, r$')
            if ti==0: ax[ci,ti].set(ylabel=r'$m \, t$')
            ax[ci,ti].tick_params(direction='in', which='both', top=True, right=True)
            ax[ci,ti].grid(ls=':', color='darkgray', alpha=0.5)

        cbar = fig.colorbar(im, ax=ax[:, :], shrink=0.6, \
                            ticks=mticker.MultipleLocator(np.pi/2), \
                            format=mticker.FuncFormatter(multiple_formatter()))
        cbar.ax.set_title(titl[0])
        plt.savefig('./plots/examples_before'+str(temp)+'.pdf', bbox_inches='tight')
        plt.show()

        fig, ax = plt.subplots(2,3, figsize = (5.6,3.4), sharey=True, sharex=True, constrained_layout=True)
        for si, sim in enumerate(dats):
            ci,ti = np.divmod(si,3)

            path2RESTsim = rest_sim_location(*exp_params, sim)
            sim, fullreal, totbeta = np.load(path2RESTsim)

            bubble2measure = fullreal[0]
            nT, nN = np.shape(bubble2measure)
            tcen, xcen = find_nucleation_center(bubble2measure, phieq, crit_thresh, crit_rad)
            tcen -= 100
            tl,tr = max(0, tcen-win), min(nT-1, tcen+win)
            xl,xr = max(0, xcen-win), min(nN-1, xcen+win)
            ext = np.array([xl-xcen,xr-xcen,tl-tcen-100,tr-tcen-100])*np.sqrt(m2(lamb))*dx

            bubble2plot = bubble2measure[tl:tr,xl:xr]
            im = ax[ci,ti].imshow(bubble2plot, interpolation=None, extent=ext, origin='lower', cmap='RdBu')

            if False:
                lavs = 6
                nT, nN = np.shape(bubble2plot)
                tt = np.linspace(tl-tcen, tr-tcen, nT)*np.sqrt(m2(lamb))*dx
                xx = np.linspace(xl-xcen, xr-xcen, nN)*np.sqrt(m2(lamb))*dx
                ttt1, xxx1 = np.meshgrid(tt, xx)
                ax[ci,ti].contour(xxx1, ttt1, bubble2plot.T, levels=lavs, aspect='auto', interpolation=None, extent=ext, origin='lower', colors='k',linewidths=0.5)

            if ci==1: ax[ci,ti].set(xlabel=r'$m \, r$')
            if ti==0: ax[ci,ti].set(ylabel=r'$m \, t$')
            ax[ci,ti].tick_params(direction='in', which='both', top=True, right=True)
            ax[ci,ti].grid(ls=':', color='darkgray', alpha=0.5)
            ax[ci,ti].legend(title=r'$v_{\rm COM}=$'+r'${:.2f}$'.format(totbeta), loc=4, fontsize=7, fancybox=True, frameon=True, framealpha=0.9, borderpad=0.15)

        cbar = fig.colorbar(im, ax=ax[:, :], shrink=0.6, \
                            ticks=mticker.MultipleLocator(np.pi/2), \
                            format=mticker.FuncFormatter(multiple_formatter()))
        cbar.ax.set_title(titl[0])
       # plt.tight_layout()
        plt.savefig('./plots/examples_after'+str(temp)+'.pdf', bbox_inches='tight')
        plt.show()

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
In [18]:
if True:
    data_photo1 = np.load('./plots/bubble_data_original_2.npy')
    rr1, ll1, rrwallfit1, llwallfit1, ttwallfit1, simulation1, ext1 = data_photo1
    data_photo2 = np.load('./plots/bubble_data_deboosted_2.npy')
    rr2, ll2, rrwallfit2, llwallfit2, ttwallfit2, simulation2, ext2 = data_photo2
    data_photo3 = np.load('./plots/bubble_data_original_1.npy')
    rr3, ll3, rrwallfit3, llwallfit3, ttwallfit3, simulation3, ext3 = data_photo3
    data_photo4 = np.load('./plots/bubble_data_deboosted_1.npy')
    rr4, ll4, rrwallfit4, llwallfit4, ttwallfit4, simulation4, ext4 = data_photo4

    rrs = np.array([rr1, rr2, rr3, rr4])
    lls = np.array([ll1, ll2, ll3, ll4])
    rrwallfits = np.array([rrwallfit1, rrwallfit2, rrwallfit3, rrwallfit4])
    llwallfits = np.array([llwallfit1, llwallfit2, llwallfit3, llwallfit4])
    ttwallfits = np.array([ttwallfit1, ttwallfit2, ttwallfit3, ttwallfit4])
    simulations = np.array([simulation1, simulation2, simulation3, simulation4])
    exts = np.array([ext1, ext2, ext3, ext4])

    exts = exts/dx2plot
    exts[:,2:]+=20
    exts = exts*dx2plot

    ttwallfits+=20

    dx2plot = dx * np.sqrt(m2(lamb))
    
    fig, ax = plt.subplots(2, 2, figsize = (5.5, 5))

    for ai, aa in enumerate(ax.flatten()):
        rr, ll = rrs[ai], lls[ai]
        rrwallfit, llwallfit, ttwallfit = rrwallfits[ai], llwallfits[ai], ttwallfits[ai]
        simulation, ext = simulations[ai], exts[ai]

        im0 = aa.imshow(simulation, interpolation='antialiased', extent=ext, origin='lower', cmap='RdBu', aspect='auto')

        aa.plot(rr*dx2plot, ttwallfit*dx2plot, color='k', ls='--', linewidth=0.5)
        aa.plot(ll*dx2plot, ttwallfit*dx2plot, color='k', ls='--', linewidth=0.5)

        aa.plot(rrwallfit*dx2plot, ttwallfit*dx2plot, color='k', ls='-', linewidth=0.7)
        aa.plot(llwallfit*dx2plot, ttwallfit*dx2plot, color='k', ls='-', linewidth=0.7)

        cbar = plt.colorbar(im0, ax=aa, shrink=0.9, \
                            ticks=mticker.MultipleLocator(np.pi/2), \
                            format=mticker.FuncFormatter(multiple_formatter()))
        cbar.ax.set_title(r'$\bar{\varphi}$')
        aa.tick_params(direction='in', which='both', top=True, right=True)
        aa.grid(ls=':', color='darkgray', alpha=0.5)

        aa.set_ylabel(r'$m \, t$', labelpad = 1)
        aa.set_xlabel(r'$m \, r$', labelpad = 1)
        aa.set_aspect(1.5)
    plt.tight_layout()
    plt.savefig('./plots/examples_deboost.pdf', dpi=500)
    plt.show()

/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:11: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  # This is added back by InteractiveShellApp.init_path()
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:12: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  if sys.path[0] == "":
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:13: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  del sys.path[0]
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:14: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:15: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  from ipykernel import kernelapp as app
In [19]:
tp = 0# 0 for average, 1 for error

titl = [r'$\bar{\varphi}$',  r'$\bar{\Pi}$', r'$\partial_r\bar{\varphi}$']

for ii, average_bubble in enumerate(bubbleList):
    fig, ax = plt.subplots(1, 3, figsize = (8,2.5))
    for cp in range(3): # 0 - field, 1 - momentum, 2 - gradient
        tmp = tmpList[ii]
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        crit_rad = 20
        crit_thresh = right_Vmax.x + 2.*sigmafld
        win = 120

        bubble2measure = average_bubble[0,0]
        nT, nN = np.shape(bubble2measure)
        tcen, xcen = find_nucleation_center(bubble2measure, phieq, crit_thresh, crit_rad)
        xcen+=1
        tl,tr = max(0, tcen-win), min(nT-1, tcen+win)
        xl,xr = max(0, xcen-win), min(nN-1, xcen+win)
        ext = np.array([xl-xcen,xr-xcen,tl-tcen,tr-tcen])*np.sqrt(m2(lamb))*dx

        bubble2plot = average_bubble[tp,cp][tl:tr,xl:xr]

        im = ax[cp].imshow(bubble2plot, interpolation=None, extent=ext, origin='lower', cmap='RdBu')
        if cp==0: cbar = plt.colorbar(im, ax=ax[cp], shrink=0.8, \
                        ticks=mticker.MultipleLocator(np.pi/2), \
                        format=mticker.FuncFormatter(multiple_formatter()))
        else:
            cbar = plt.colorbar(im, ax=ax[cp], shrink=0.8)
            lls = cbar.ax.get_ylim()
            a = cbar.ax.get_yticks().tolist()
            a = [round(al,2) for ja,al in enumerate(a) if ja%2==1][::]
            cbar.ax.set_yticks(a)
            #cbar.ax.set_ylim(a[0], a[-1])
            cbar.ax.set_ylim(lls)
            a = [r'${:.1f}$'.format(al) for al in a]
            cbar.ax.set_yticklabels(a)
            
        cbar.ax.set_title(titl[cp])

        nT, nN = np.shape(bubble2plot)
        tt = np.linspace(tl-tcen, tr-tcen, nT)*np.sqrt(m2(lamb))*dx
        xx = np.linspace(xl-xcen, xr-xcen, nN)*np.sqrt(m2(lamb))*dx
        ttt1, xxx1 = np.meshgrid(tt, xx)

        lavs = [8, 6, 8][cp]
        ax[cp].grid(True, color='k', linewidth=0.5, ls=':', alpha=0.2)
        ax[cp].contour(xxx1, ttt1, bubble2plot.T, levels=lavs, aspect='auto', interpolation=None, extent=ext, origin='lower', colors='k',linewidths=0.5)

        labss = r'$T/m = {}$'.format(round(temp/np.sqrt(m2(lamb)),1))
        
        ax[cp].grid(ls=':', color='darkgray', alpha=0.7)
        ax[cp].xaxis.set_minor_locator(MultipleLocator(1))
        ax[cp].yaxis.set_minor_locator(MultipleLocator(1))
        ax[cp].set_xlabel(r'$m \, r$', labelpad = 1)
        ax[0].set_ylabel(r'$m \, t$', labelpad = 1)

        ax[cp].yaxis.set_ticks_position('both')
        ax[cp].xaxis.set_ticks_position('both')
        ax[cp].tick_params(which='both', axis="y", direction="in")
        ax[cp].tick_params(which='both', axis="x", direction="in")
    plt.tight_layout()
    plt.savefig('./plots/average_bubble'+str(temp)+'.pdf', bbox_inches='tight')
    plt.show()

/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:50: UserWarning: The following kwargs were not used by contour: 'aspect', 'interpolation'
In [20]:
if False:
    for ii, average_bubble in enumerate(bubbleList):
        average_bubble = bubbleList[ii]
        tmp = tmpList[ii]
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        print(*exp_params)

        critical_sim = extract_data(nLat, critical_sim_file(*exp_params))
        precursor_sim = extract_data(nLat, precursor_sim_file(*exp_params))

        titls = [r'$\rm Critical$', r'$\rm Precursor$']

        tdecay = max(0, np.shape(critical_sim)[1] - nLat//2)
        crit_thresh = right_Vmax.x + 5.*sigmafld
        critical_sim, crit_rad = centre_bubble(critical_sim, tdecay, phieq, crit_thresh)
        critical_sim = remove_collisions(critical_sim, phieq)

        print('Decay time standard, radius:', tdecay, crit_rad)
        critical_sim = critical_sim[0]
        precursor_sim = precursor_sim[0]

     #   critical_sim = gaussian_filter(critical_sim, sigma=5, mode='nearest')
     #   precursor_sim = gaussian_filter(precursor_sim, sigma=5, mode='nearest')

        truebubble = critical_sim[:nLat//4, 3*nLat//8 : 5*nLat//8]
        precursor  = precursor_sim[:nLat//4, 3*nLat//8 : 5*nLat//8]
        print(np.shape(precursor))


        nN, nT      = np.shape(truebubble)
        tcen, xcen  = find_nucleation_center(truebubble, phieq, crit_thresh, crit_rad//3)
        xcen += 1
        tl,tr,xl,xr = -tcen, nT-1-tcen, -xcen, nN-1-xcen
        exts1       = np.array([xl, xr, tl, tr])*dx*np.sqrt(m2(lamb))
        tt1, xx1    = np.linspace(tl, tr, nT)*dx*np.sqrt(m2(lamb)), np.linspace(xl, xr, nN)*dx*np.sqrt(m2(lamb))
        ttt1, xxx1  = np.meshgrid(tt1, xx1)

        tl,tr,xl,xr = 0, nT-1, -xcen, nN-1-xcen
        exts2       = np.array([xl, xr, tl, tr])*dx*np.sqrt(m2(lamb))
        tt2, xx2    = np.linspace(tl, tr, nT)*dx*np.sqrt(m2(lamb)), np.linspace(xl, xr, nN)*dx*np.sqrt(m2(lamb))
        ttt2, xxx2  = np.meshgrid(tt2, xx2)

        fig, ax = plt.subplots(1, 2, figsize = (6,3))
        ax[0].contour(xxx1, ttt1, truebubble.T, levels=3, aspect='auto', interpolation=None, extent=exts1, origin='lower', colors='k', linewidths=0.5)
        im0 = ax[0].imshow(truebubble, aspect='auto', interpolation='none', extent=exts1, origin='lower', cmap='RdBu')
        clb0 = plt.colorbar(im0, ax = ax[0], ticks = mticker.MultipleLocator(np.pi/2), format = mticker.FuncFormatter(multiple_formatter()))
        clb0.ax.set_title(r'$\bar{\varphi}$')

        ax[1].contour(xxx2, ttt2, np.abs(np.pi-precursor.T), levels=3, aspect='auto', interpolation=None, extent=exts2, origin='lower', colors='k', linewidths=0.5)
        im1 = ax[1].imshow(precursor, aspect='auto', interpolation='none', extent=exts2, origin='lower', cmap='RdBu')
        clb1 = plt.colorbar(im1, ax = ax[1], ticks = mticker.MultipleLocator(np.pi/2), format = mticker.FuncFormatter(multiple_formatter()))
        clb1.ax.set_title(r'$\bar{\varphi}$')

        for ai, aa in enumerate(ax):
          #  aa.text(5.5,.5,'TEST TEST TEST TEST',
          #            bbox={'facecolor':'white','alpha':1,'edgecolor':'none','pad':1},
          #            ha='center', va='center') 
            aa.text([6.8,6.2][ai], [12.,18.][ai], titls[ai], ha='center', va='center', \
                    bbox={'boxstyle':'round','facecolor':'white','alpha':0.85,'edgecolor':'k','pad':0.3}, fontsize=10)

            aa.grid(ls=':', color='darkgray', alpha=0.7)
            aa.xaxis.set_minor_locator(MultipleLocator(1))
            aa.yaxis.set_minor_locator(MultipleLocator(1))
          #  aa.set_xlabel(r'$\phi_0^{-1} \sqrt{V_0} \; r$')
          #  aa.set_ylabel(r'$\phi_0^{-1} \sqrt{V_0} \; t$')
            aa.set_xlabel(r'$m \, r$')
            aa.set_ylabel(r'$m \, t$')
            aa.yaxis.set_ticks_position('both')
            aa.xaxis.set_ticks_position('both')
            aa.tick_params(which='both', axis="y", direction="in")
            aa.tick_params(which='both', axis="x", direction="in")
           # a = aa.get_xticks().tolist()[1:-1:]
           # a = [round(al,2) for al in a]
           # aa.set_xticks(a)
           # aa.set_xticklabels(a)
            a = aa.get_yticks().tolist()[1:-1:2]
            a = [round(al,2) for al in a]
            aa.set_yticks(a)
            a = [r'${:.0f}$'.format(al) for al in a]
            aa.set_yticklabels(a)

        plt.tight_layout()
        plt.savefig('./plots/critical_and_precursor'+str(temp)+'.pdf', bbox_inches='tight', dpi=500)
        plt.show()
In [21]:
if False:
    titls = [r'$\rm Critical$', r'$\rm Subcritical$']

    super_instanton_sim = extract_data(nLat, supercrit_instanton_file(*exp_params))
    instanton_sim = extract_data(nLat, instanton_file(*exp_params))

    tdecay = max(0, np.shape(super_instanton_sim)[1] - nLat//2)
    crit_thresh = right_Vmax.x + 5.*sigmafld
    super_instanton_sim, crit_rad = centre_bubble(super_instanton_sim, tdecay, phieq, crit_thresh)
    super_instanton_sim = remove_collisions(super_instanton_sim, phieq)

    for tmp in range(len(tempList)):
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

     #   np.save(supercrit_instanton_file(*exp_params), super_instanton_sim[0,0,:])
     #   np.save(instanton_file(*exp_params), instanton_sim[0,0,:])

    print('Decay time standard, radius:', tdecay, crit_rad)
    super_bubble = super_instanton_sim[0, 40:nLat//3+41, nLat//3:2*nLat//3+1]
    bubble       = instanton_sim[0, 40:nLat//3+41, nLat//3:2*nLat//3+1]

    nN, nT      = np.shape(super_bubble)
    tcen, xcen  = find_nucleation_center(super_bubble, phieq, crit_thresh, crit_rad//3)

    tl,tr,xl,xr = -tcen, nT-1-tcen, -xcen, nN-1-xcen
    exts1       = np.array([xl, xr, tl, tr])*dx*np.sqrt(m2(lamb))
    tt1, xx1    = np.linspace(tl, tr, nT)*dx*np.sqrt(m2(lamb)), np.linspace(xl, xr, nN)*dx*np.sqrt(m2(lamb))
    ttt1, xxx1  = np.meshgrid(tt1, xx1)

    tl,tr,xl,xr = 0, nT-1, -xcen, nN-1-xcen
    exts2       = np.array([xl, xr, tl, tr])*dx*np.sqrt(m2(lamb))
    tt2, xx2    = np.linspace(tl, tr, nT)*dx*np.sqrt(m2(lamb)), np.linspace(xl, xr, nN)*dx*np.sqrt(m2(lamb))
    ttt2, xxx2  = np.meshgrid(tt2, xx2)

    fig, ax = plt.subplots(1, 2, figsize = (6,3))
    ax[0].contour(xxx1, ttt1, super_bubble.T, levels=4, aspect='auto', interpolation='none', extent=exts1, origin='lower', colors='k', linewidths=0.5)
    im0 = ax[0].imshow(super_bubble, aspect='auto', interpolation='none', extent=exts1, origin='lower', cmap='RdBu')
    clb0 = plt.colorbar(im0, ax = ax[0], ticks = mticker.MultipleLocator(np.pi/2), format = mticker.FuncFormatter(multiple_formatter()))
    clb0.ax.set_title(r'$\bar{\varphi}$')

    ax[1].contour(xxx2, ttt2, bubble.T, levels=7, aspect='auto', interpolation='none', extent=exts2, origin='lower', colors='k', linewidths=0.5)
    im1 = ax[1].imshow(bubble, aspect='auto', interpolation='none', extent=exts2, origin='lower', cmap='RdBu')
    clb1 = plt.colorbar(im1, ax = ax[1], ticks = mticker.MultipleLocator(np.pi/2), format = mticker.FuncFormatter(multiple_formatter()))
    clb1.ax.set_title(r'$\bar{\varphi}$')

    for ai, aa in enumerate(ax):
      #  aa.text(5.5,.5,'TEST TEST TEST TEST',
      #            bbox={'facecolor':'white','alpha':1,'edgecolor':'none','pad':1},
      #            ha='center', va='center') 
        aa.text([9.,7.2][ai], [-7.8,2.][ai], titls[ai], ha='center', va='center', \
                bbox={'boxstyle':'round','facecolor':'white','alpha':0.85,'edgecolor':'k','pad':0.3}, fontsize=9)

        aa.grid(ls=':', color='darkgray', alpha=0.7)
        aa.xaxis.set_minor_locator(MultipleLocator(1))
        aa.yaxis.set_minor_locator(MultipleLocator(1))
      #  aa.set_xlabel(r'$\phi_0^{-1} \sqrt{V_0} \; r$')
      #  aa.set_ylabel(r'$\phi_0^{-1} \sqrt{V_0} \; t$')
        aa.set_xlabel(r'$m \, r$')
        aa.set_ylabel(r'$m \, t$')
        aa.yaxis.set_ticks_position('both')
        aa.xaxis.set_ticks_position('both')
        aa.tick_params(which='both', axis="y", direction="in")
        aa.tick_params(which='both', axis="x", direction="in")
    plt.tight_layout()
    plt.savefig('./plots/bare_spheleron.pdf', dpi=500)
    plt.show()
    
    tmp = 0
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])
    print(*exp_params)

    instanton = np.load(instanton_file(*exp_params))# - np.pi
    supercritinstanton = np.load(supercrit_instanton_file(*exp_params))# - np.pi
    print(np.shape(instanton), np.shape(supercritinstanton))

    xmin, xmax = nLat//3, nLat*2//3
    #plt.plot(xlist[xmin:xmax], instanton[xmin:xmax])
    #plt.plot(xlist[xmin:xmax], supercritinstanton[xmin:xmax])
    plt.plot(xlist[xmin:xmax], supercritinstanton[xmin+1:xmax+1] - instanton[xmin:xmax])
    plt.plot(xlist[xmin:xmax], supercritinstanton[xmin:xmax] - super_instanton_sim[0,0,xmin:xmax])
    plt.plot(xlist[xmin:xmax], instanton[xmin:xmax] - instanton_sim[0,0,xmin:xmax])
    plt.show()

    crit_soln_f90 = '(/'
    for ind, iii in enumerate(instanton):
        crit_soln_f90 = crit_soln_f90 + str(iii)
        if ind != len(instanton)-1:
            crit_soln_f90 = crit_soln_f90 + ', '
    crit_soln_f90 += '/)'
    #print(crit_soln_f90)
In [22]:
if False:
    tmp = 5
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])

    undecayed_sims = np.load(sims_notdecayed_file(*exp_params, minSim, maxSim, nTimeMAX))
    print(undecayed_sims[:,0])

    decay_times  = np.load(decay_times_file(*exp_params, minSim, maxSim, nTimeMAX))
    
    minDecTime  = 4*nLat
    alltimes    = decay_times[:,1]
    simList2Do  = decay_times[alltimes>=minDecTime, 0]
    print(simList2Do)

  #  sims_decayed   = np.load(sims_decayed_file(*exp_params, minSim, maxSim, nTimeMAX))

    choose = random.sample(undecayed_sims[:,0].tolist(), 1)
    for sim in choose:
        path2sim = sim_location(*exp_params, sim)
        
        real, cho = get_realisation(nLat, sim, phieq, path2sim)
        nC, nT, nN = np.shape(real)
        print('Sim, nT, outcome', sim, nT, cho)

        bubble = real[0,:]#nT-3*nLat//4]
        nT, nN = np.shape(bubble)
        ax1 = simple_imshow(bubble, [0,nT,0,nN], 'Simulation '+str(sim))
        plt.show()
        
        fftbubble = np.fft.fft2(bubble)
        nT, nN = np.shape(fftbubble)
        #fftbubble = fftbubble[:nT//2,:nN//2]
        fftbubble[:, nN*4//10:nN*6//10+1] = 0.
        ax2 = simple_imshow(np.log(fftbubble.real), [0,nT,0,nN], 'Simulation '+str(sim))

        bubble2 = np.fft.ifft2(fftbubble)
        nT, nN = np.shape(bubble2)
        ax3 = simple_imshow(bubble2.real, [0,nT,0,nN], 'Simulation '+str(sim))
        ax4 = simple_imshow(bubble2.imag, [0,nT,0,nN], 'Simulation '+str(sim))

        ax5 = simple_imshow(bubble2.real-bubble, [0,nT,0,nN], 'Simulation '+str(sim))
        
        if False:
#        if sim in sims_decayed[:,0]:
            path2CLEANsim = clean_sim_location(*exp_params, sim)
            real, sim, tdecay, outcome = np.load(path2CLEANsim)

            bubble = real[0,-nLat:]
            nT, nN = np.shape(bubble)
            simple_imshow(bubble, [0,nT,0,nN], 'Simulation '+str(sim))
            print('Sim, outcome0, outcome1, tdecay', sim, outcome, cho, tdecay)
In [23]:
# Classify decays
fig, ax = plt.subplots(1, len(tempList), figsize = (15,2.5))
plt.style.use('seaborn-whitegrid') # nice and clean grid
cols = cycle(['purple', 'forestgreen', 'orange', 'blue'])

for tmp in range(len(tempList)):
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])

    undecayed_sims = np.load(sims_notdecayed_file(*exp_params, minSim, maxSim, nTimeMAX))
    decayed_sims = np.load(sims_decayed_file(*exp_params, minSim, maxSim, nTimeMAX))
    decay_times = np.load(decay_times_file(*exp_params, minSim, maxSim, nTimeMAX))

    outcomes = decayed_sims[:,1]; print(len(outcomes))

    labs = r'$T/m={:.1f}$'.format(temp/np.sqrt(m2(lamb)))
    ax[tmp].hist(outcomes, density=True, bins=2, rwidth=0.8, color=next(cols), alpha=0.9, linewidth=0.5)
    ax[tmp].set_xlabel(r'$2\pi \quad\quad 0$')
    ax[tmp].set_ylabel('PDF')
    ax[tmp].set_title(labs)
plt.tight_layout()
plt.savefig('./plots/vacuum_choice.pdf', rasterize=True); plt.show()

1003
2192
3165
3711
3910
3984

/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:22: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "rasterize" which is no longer supported as of 3.3 and will become an error in 3.6
In [24]:
tmp = tmpList[0]
phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)

instanton = instantonList[tmp]
gradinst  = (np.roll(instanton,-1) - instanton)/dx2plot

testgrd = fldcritList[tmp][-1,:]*dx/dx2plot

testfld = fldcritList[tmp][0,:]
testfld = (np.roll(testfld,-1) - testfld)/dx2plot

xtest = (np.arange(nLat) - nLat//2) * dx2plot
cds   = np.array(np.linspace(nLat//2-len(testgrd)//2, nLat//2+len(testgrd)//2, len(testgrd)), dtype='int')

fig, ax = plt.subplots(1, 1, figsize = (5,3.5))
ax.plot(xtest[cds], testgrd, ls='-', label='gradient of average bubble')
ax.plot(xtest[cds], testfld, ls='-', label='average gradient of average bubble')
ax.plot(xtest[cds], gradinst[cds], ls='-',color='k', label='gradient of critical lattice solution')

ax.legend()
ax.set_ylabel(r'$\partial_x \phi_{\rm crit}(x)$')
ax.set_xlabel(r'$\phi_0^{-1} \sqrt{V_0} \; x$')
plt.grid(True, ls=':', color='lightgray', alpha=0.5)
plt.show()
In [25]:
if True:
    for tmp in range(len(tempList)):
        if tmp==4: break
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        # 0 for field; 1 for momentum 
        find = 0
        # modes to plot:
        aa, bb = 1, knyq-1
        # time steps to plot
        tstep = 10
        # if tind not in np.array([0,2,7,27,46]): continue

        ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
        tlist, PSfld1 = ALL_powspec1[0][::tstep], ALL_powspec1[1][:, find, ::tstep, aa:bb]
        del ALL_powspec1

        ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
        tlist, PSfld2 = ALL_powspec2[0][::tstep], ALL_powspec2[1][:, find, ::tstep, aa:bb]
        del ALL_powspec2

        fig, ax = plt.subplots(1,1, figsize = (6,4.))
        plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), pspec(lamb, phi0, temp)[aa:bb], color='k', linewidth=0.5)
        for tind, tt in enumerate(tlist):
            curve = np.nanmean(np.concatenate((PSfld1[:,tind], PSfld2[:,tind]), axis=0), axis=0)

            lab = r'${:.1f}$'.format(round(tt*dx*np.sqrt(m2(lamb)),1))
            plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), curve, ls='-', label=lab)

        del PSfld1, PSfld2

        h, l = ax.get_legend_handles_labels() # Extracting handles and labels
        handles = [plt.plot([],marker="", ls="")[0]] + h
        labels = [r'$m\,t:$'] + l  # Merging labels
        leg = ax.legend(handles, labels, ncol=2, frameon=False, loc=1)

        ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
        ax.tick_params(direction='in', which='both', top=True, right=True)
        ax.set_xscale('log')
        ax.set_ylabel(r'$\left\langle \left|\bar{\varphi}_k(t)\right|^2 \right\rangle $')
        ax.set_xlabel(r'$k/m$')
        plt.tight_layout()
        plt.savefig('./plots/powespec_tevol.pdf')
        plt.show()
In [26]:
if True:
    for tmp in range(len(tempList)):
        if tmp==4: break
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        # 0 for field; 1 for momentum 
        find = 1
        # modes to plot:
        aa, bb = 1, knyq-1
        # time steps to plot
        tstep = 10

        ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
        tlist, PSfld1 = ALL_powspec1[0][::tstep], ALL_powspec1[1][:, find, ::tstep, aa:bb]
        del ALL_powspec1

        ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
        tlist, PSfld2 = ALL_powspec2[0][::tstep], ALL_powspec2[1][:, find, ::tstep, aa:bb]
        del ALL_powspec2

        fig, ax = plt.subplots(1,1, figsize = (4.7,3.))
        plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), (w2(lamb)*pspec(lamb, phi0, temp))[aa:bb], color='k', linewidth=0.5)

        for tind, tt in enumerate(tlist):
            curve = np.nanmean(np.concatenate((PSfld1[:,tind], PSfld2[:,tind]), axis=0), axis=0)

            lab = r'${:.1f}$'.format(round(tt*dx*np.sqrt(m2(lamb)),1))
            plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), curve, ls='-', label=lab)

        del PSfld1, PSfld2

        h, l = ax.get_legend_handles_labels() # Extracting handles and labels
        handles = [plt.plot([],marker="", ls="")[0]] + h
        labels = [r'$m\,t:$'] + l  # Merging labels
        leg = ax.legend(handles, labels, ncol=3, frameon=False, bbox_to_anchor=(1.5, 0.5))

        ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
        ax.tick_params(direction='in', which='both', top=True, right=True)
        ax.set_xscale('log')
        ax.set_ylabel(r'$\left\langle \left|\bar{\Pi}_k(t)\right|^2 \right\rangle $')
        ax.set_xlabel(r'$k/m$')
        plt.tight_layout()
        plt.savefig('./plots/powespec_tevol.pdf')
        plt.show()

/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:43: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all axes decorations.
In [27]:
if True:
    savemeff = np.zeros((4))
    for tmp in range(len(tempList)):
        if tmp==4: break
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        # modes to plot:
        aa, bb = 1, knyq//5

        ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
        tlist, PSfld1, PSmom1 = ALL_powspec1[0], ALL_powspec1[1][:, 0, -50:, aa:bb], ALL_powspec1[1][:, 1, -50:, aa:bb]
        del ALL_powspec1

        ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
        tlist, PSfld2, PSmom2 = ALL_powspec2[0], ALL_powspec2[1][:, 0, -50:, aa:bb], ALL_powspec2[1][:, 1, -50:, aa:bb]
        del ALL_powspec2

        avPSfld = np.nanmean(np.nanmean(np.concatenate((PSfld1, PSfld2), axis=0), axis=1), axis=0)
        avPSmom = np.nanmean(np.nanmean(np.concatenate((PSmom1, PSmom2), axis=0), axis=1), axis=0)
        del PSfld1, PSfld2, PSmom1, PSmom2
        curve   = (avPSmom/avPSfld)

        norml    = 1./ phi0 / np.sqrt(2. * lenLat)
        w2m      = lambda ks, m: ks**2. + m**2.
        pofk     = lambda ks, m: w2m(ks, m)
        pred_fit = lambda x, data: sco.curve_fit(pofk, x, data, p0=0.09)[0]

        fig, ax = plt.subplots(1,1, figsize = (4.7,3.))
        plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), pofk(klist[aa:bb], np.sqrt(m2(lamb))), color='k', lw=0.5, ls='-')

        lab1 = r'$m={:.4f}$'.format(np.sqrt(m2(lamb)), temp)
        plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), curve, ls='-', label=lab1)

        best_ps = pred_fit(klist[aa:bb], curve)
        lab2 = r'$m={:.4f}$'.format(*best_ps)

        savemeff[tmp] = best_ps[0]
        plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), pofk(klist[aa:bb], *best_ps), color='r', lw=1, ls=':', label=lab2)

        leg = ax.legend(ncol=1, frameon=False, loc='best')#, bbox_to_anchor=(2.3, 1.))
        ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
        ax.tick_params(direction='in', which='both', top=True, right=True)
        ax.set_xscale('log')
        ax.set_yscale('log')
        ax.set_ylabel(r'$\left\langle \left|\bar{\Pi}_k(t)\right|^2/\left|\bar{\varphi}_k(t)\right|^2 \right\rangle $')
        ax.set_xlabel(r'$k/m$')
        plt.tight_layout()
        #plt.savefig('./plots/bestfit_masses.pdf')
        plt.show()
        
print(savemeff)

/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:19: RuntimeWarning: Mean of empty slice
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:20: RuntimeWarning: Mean of empty slice

[0.08631948 0.08449074 0.08163867 0.08173309]
In [28]:
if True:
    norml    = 1./ phi0 / np.sqrt(2. * lenLat)
    w2m      = lambda ks, msq: ks**2. + msq
    pofk0    = lambda ks, msq, te: norml / w2m(ks, msq)**0.25 * np.sqrt(2./(np.exp(w2m(ks, msq)**0.5/te) - 1.))
    #pofk     = lambda ks, msq, te: norml / w2m(ks, msq)**0.25 * np.sqrt(2./(w2m(ks, msq)**0.5/te + 0.5*(w2m(ks, msq)**0.5/te)**2. ))
    pofk     = lambda ks, msq, te: norml / w2m(ks, msq)**0.25 * np.sqrt(2./(w2m(ks, msq)**0.5/te))
    f_pred   = lambda ks, msq, te: np.abs(pofk(ks, msq, te))**2.
    f_pred1  = lambda ks, te: np.abs(pofk(ks, savemeff[tmp]**2., te))**2.
    f_pred2  = lambda ks, te: np.abs(pofk0(ks, m2(lamb), te))**2.
    pred_fit = lambda x, data: sco.curve_fit(f_pred1, x, data)[0]

    save_teff_from_fld = np.zeros((4))
    for tmp in range(len(tempList)):
        if tmp==4: break
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        # modes to plot:
        aa, bb = 1, knyq//5
        # time slices to average
        am, bm = -50, -1

        ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
        tlist, PSfld1 = ALL_powspec1[0], ALL_powspec1[1][:, 0, am:bm, aa:bb]
        del ALL_powspec1

        ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
        tlist, PSfld2 = ALL_powspec2[0], ALL_powspec2[1][:, 0, am:bm, aa:bb]
        del ALL_powspec2

        PSfld = np.concatenate((PSfld1, PSfld2), axis=0)
        curve = np.nanmean(np.nanmean(PSfld, axis=1), axis=0)
        del PSfld1, PSfld2

        fig, ax = plt.subplots(1,1, figsize = (4.7,3.))
        plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), f_pred2(klist[aa:bb], temp), color='k', lw=1, ls='-', label=r'$\rm BE, \; initially$')

        lab1 = r'${:.1f} < m\,t $'.format(round(tlist[am]*dx*np.sqrt(m2(lamb)),1)) + r'$< {:.1f}$'.format(round(tlist[bm]*dx*np.sqrt(m2(lamb)),1))
        [plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), ii, ls='-', label=(lab1 if i0==0 else None)) for i0, ii in enumerate(np.nanmean(PSfld, axis=0))]
        #plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), curve, ls='-', lw=2, label=lab1)

        best_ps = pred_fit(klist[aa:bb], curve)
        lab2 = r'$T={:.4f}$'.format(*best_ps)
        save_teff_from_fld[tmp] = best_ps[-1]
        plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), f_pred1(klist[aa:bb], *best_ps), color='r', lw=1, ls=':', label=lab2)

        leg = ax.legend(title=r'$T={:.2f}$'.format(temp), ncol=1, frameon=False, loc='best')
        ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
        ax.tick_params(direction='in', which='both', top=True, right=True)
        ax.set_xscale('log')
        ax.set_ylabel(r'$\left\langle \left|\bar{\varphi}_k(t)\right|^2 \right\rangle $')
        ax.set_xlabel(r'$k/m$')
        plt.tight_layout()
        #plt.savefig('./plots/powespec_tevol.pdf')
        plt.show()

/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:32: RuntimeWarning: Mean of empty slice
In [29]:
if True:
    norml    = 1./ phi0 / np.sqrt(2. * lenLat)
    w2m      = lambda ks, msq: ks**2. + msq
    pofk0    = lambda ks, msq, te: norml * w2m(ks, msq)**0.25 * np.sqrt(2./(np.exp(w2m(ks, msq)**0.5/te) - 1.))
    #pofk     = lambda ks, msq, te: norml * w2m(ks, msq)**0.25 * np.sqrt(2./(w2m(ks, msq)**0.5/te + 0.5*(w2m(ks, msq)**0.5/te)**2. ))
    pofk     = lambda ks, msq, te: norml * w2m(ks, msq)**0.25 * np.sqrt(2./(w2m(ks, msq)**0.5/te))
    f_pred   = lambda ks, msq, te: np.abs(pofk(ks, msq, te))**2.
    f_pred1  = lambda ks, te: np.abs(pofk(ks, savemeff[tmp]**2., te))**2.
    f_pred2  = lambda ks, te: np.abs(pofk0(ks, m2(lamb), te))**2.
    pred_fit = lambda x, data: sco.curve_fit(f_pred1, x, data)[0]

    save_teff_from_mom = np.zeros((4))
    for tmp in range(len(tempList)):
        if tmp==4: break
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        # modes to plot:
        aa, bb = 1, knyq//5
        # time slices to average
        am, bm = -50, -1

        ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
        tlist, PSfld1 = ALL_powspec1[0], ALL_powspec1[1][:, 1, am:bm, aa:bb]
        del ALL_powspec1

        ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
        tlist, PSfld2 = ALL_powspec2[0], ALL_powspec2[1][:, 1, am:bm, aa:bb]
        del ALL_powspec2

        PSfld = np.concatenate((PSfld1, PSfld2), axis=0)
        curve = np.nanmean(np.nanmean(PSfld, axis=1), axis=0)
        del PSfld1, PSfld2

        fig, ax = plt.subplots(1,1, figsize = (4.7,3.))
        #plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), f_pred2(klist[aa:bb], temp), color='k', lw=1, ls='-', label=r'$\rm BE, \; initially$')

        lab1 = r'${:.1f} < m\,t $'.format(round(tlist[am]*dx*np.sqrt(m2(lamb)),1)) + r'$< {:.1f}$'.format(round(tlist[bm]*dx*np.sqrt(m2(lamb)),1))
        [plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), ii, ls='-', label=(lab1 if i0==0 else None)) for i0, ii in enumerate(np.nanmean(PSfld, axis=0))]
        #plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), curve, ls='-', lw=2, label=lab1)

        best_ps = pred_fit(klist[aa:bb//10], curve[:bb//10-1])
        lab2 = r'$T={:.4f}$'.format(*best_ps)
        save_teff_from_mom[tmp] = best_ps[-1]
        plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), f_pred1(klist[aa:bb], *best_ps), color='r', lw=1, ls=':', label=lab2)

        lab3 = r'$T={:.4f}$'.format(save_teff_from_fld[tmp])
        plt.plot(klist[aa:bb]/np.sqrt(m2(lamb)), f_pred1(klist[aa:bb], save_teff_from_fld[tmp]), color='g', lw=1, ls=':', label=lab3)

        leg = ax.legend(title=r'$T={:.2f}$'.format(temp), ncol=1, frameon=False, loc='best')
        ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
        ax.tick_params(direction='in', which='both', top=True, right=True)
        ax.set_xscale('log')
        ax.set_ylabel(r'$\left\langle \left|\bar{\Pi}_k(t)\right|^2 \right\rangle $')
        ax.set_xlabel(r'$k/m$')
        plt.tight_layout()
        #plt.savefig('./plots/powespec_tevol.pdf')
        plt.show()

/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:32: RuntimeWarning: Mean of empty slice
In [30]:
fig, ax = plt.subplots(1,1, figsize = (4.7,3.))
plt.plot(tempList[tmpList], save_teff_from_fld, label=r'$\rm Best \; Fit \; RJ \; from Field$')
plt.plot(tempList[tmpList], save_teff_from_mom, label=r'$\rm Best \; Fit \; RJ \; from Momentum$')
plt.plot(tempList[tmpList], tempList[tmpList], label=r'$\rm Bare \; Temperature$')
plt.legend()
plt.show()
In [31]:
exp_bounce = lambda x, a, b: a * np.sqrt(b/x) * np.exp(-b/x)
get_best_bounce = lambda x, y: sco.curve_fit(exp_bounce, x, y)[0]
In [32]:
fig, ax = plt.subplots(1, 3, figsize = (9.7,3.), gridspec_kw={'width_ratios': [4, 2.3, 2.3]})
gammas1 = np.zeros((len(tempList), 2))
cls = cycle(allcolors[:4])

for tmp in range(len(tempList[:4])):
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])
    labss = r'${}$'.format(round(temp/np.sqrt(m2(lamb)),1))

    lmmax, lmmin = 5*nLat, 1.5*nLat#2*nLat/3

    decay_times = np.load(decay_times_file(*exp_params, minSim, maxSim, nTimeMAX))
    ndcys = maxSim - minSim
    decay_times = np.sort(decay_times[:,1])
    decay_times = decay_times[decay_times < lmmax]

    ax[0].fill_betweenx(np.linspace(-1,2,10), lmmin*dx*np.sqrt(m2(lamb)), lmmax*dx*np.sqrt(m2(lamb)), color='darkgray', alpha=0.05)

    frmin = [np.argmin(np.abs(decay_times - lmmax))]
    frmax = [np.argmin(np.abs(decay_times - lmmin))]
    frmin = survive_prob(decay_times, ndcys)[frmin]
    frmax = survive_prob(decay_times, ndcys)[frmax]

    tmin = decay_times[np.argmin(np.abs(survive_prob(decay_times, ndcys) - frmax))]
    tmax = decay_times[np.argmin(np.abs(survive_prob(decay_times, ndcys) - frmin))]
    jfit_times   = lin_fit_times(decay_times, ndcys, tmin, tmax)
    gammas1[tmp] = np.array([temp, -jfit_times[0]])

    col = next(cls)
    ax[0].plot(decay_times*dx*np.sqrt(m2(lamb)), np.exp(get_line(decay_times, *jfit_times)), ls=':', lw=1, color='k')
    ax[0].plot(decay_times*dx*np.sqrt(m2(lamb)), survive_prob(decay_times, ndcys), color=col, ls='-', label=labss)

temprange = np.linspace(tempList[0], tempList[3], 50)
bounce_dat0 = get_best_bounce(tempList[:4]/np.sqrt(m2(lamb)), gammas1[:4,1]/np.sqrt(m2(lamb)))
ax[1].plot(temprange/np.sqrt(m2(lamb)), exp_bounce(temprange/np.sqrt(m2(lamb)), *bounce_dat0), \
           color='k', label=r'$A \sqrt{\frac{E}{T}} e^{-E/T}$')

temprange = np.linspace(save_teff_from_fld[0], save_teff_from_fld[3], 50)
bounce_dat = get_best_bounce(save_teff_from_fld[:4]/np.sqrt(m2(lamb)), gammas1[:4,1]/np.sqrt(m2(lamb)))
ax[2].plot(temprange/np.sqrt(m2(lamb)), exp_bounce(temprange/np.sqrt(m2(lamb)), *bounce_dat), \
           color='k', label=r'$A \sqrt{\frac{E}{T_{\rm eff}}} e^{-E/T_{\rm eff}}$')
print(bounce_dat)

cls = cycle(allcolors[:4])
for tmp in range(len(tempList)):
    if tmp>3: break
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])

    col = next(cls)
    ax[1].plot(temp/np.sqrt(m2(lamb)), gammas1[tmp,1]/np.sqrt(m2(lamb)), color=col, marker='o', linestyle='none')
    
    temp2 = save_teff_from_fld[tmp]
    ax[2].plot(temp2/np.sqrt(m2(lamb)), gammas1[tmp,1]/np.sqrt(m2(lamb)), color=col, marker='o', linestyle='none')

ax[0].set_yscale('log')
#ax[0].ticklabel_format(axis='both', scilimits=[-5,5.])
ax[0].set_xlim(0, lmmax*dx*np.sqrt(m2(lamb)))
ax[0].set_ylim(4e-2, 1.05)
ax[0].set_ylabel(r'$\rm Pr(survive)$', fontsize=13)
ax[0].set_xlabel(r'$m \, t$')

h, l = ax[0].get_legend_handles_labels() # Extracting handles and labels
handles = [plt.plot([],marker="", ls="")[0]] + h
labels = [r'$T/m:$'] + l  # Merging labels
leg = ax[0].legend(handles, labels, ncol=5, loc='lower left', labelspacing=0.3, columnspacing=0.6)#, markerfirst=False)
ax[0].grid(True, ls=':', color='lightgray', alpha=0.5)
ax[0].tick_params(direction='in', which='both', top=True, right=True)

for vpack in leg._legend_handle_box.get_children()[:1]:
    for hpack in vpack.get_children():
        for ii in hpack.get_children():
            ii.set_width(-10)
for vpack in leg._legend_handle_box.get_children()[1:]:
    for hpack in vpack.get_children():
        for ii in hpack.get_children():
            ii.set_width(18)


a = tempList[tmpList]/np.sqrt(m2(lamb))
ax[1].set_xticks(a)

a = [round(ii,2) for ii in save_teff_from_fld/np.sqrt(m2(lamb))]
ax[2].set_xticks(a)

for ai, aa in enumerate(ax[1:]):
    aa.set_ylabel(r'$\Gamma / m \, L$', fontsize=13)
    aa.set_xlabel([r'$T/m$',r'$T_{\rm eff}/m$'][ai])
    aa.legend(loc=2)
    aa.grid(True, ls=':', color='lightgray', alpha=0.5)
    aa.tick_params(direction='in', which='both', top=True, right=True)
    aa.ticklabel_format(axis='both', style='scientific', scilimits=[-1.,0.])

fig.tight_layout()
plt.savefig('./plots/decay_rate_and_surv_fraction.pdf')
plt.show()

[1.094875   3.09474548]
In [33]:
fig, ax = plt.subplots(1, 3, figsize = (9.7,3.), gridspec_kw={'width_ratios': [4, 2.3, 2.3]})
gammas1 = np.zeros((len(tempList), 2))
cls = cycle(allcolors[:4])

for tmp in range(len(tempList[:4])):
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])
    labss = r'${}$'.format(temp)

    lmmax, lmmin = 5*nLat, 1.5*nLat#2*nLat/3

    decay_times = np.load(decay_times_file(*exp_params, minSim, maxSim, nTimeMAX))
    ndcys = maxSim - minSim
    decay_times = np.sort(decay_times[:,1])
    decay_times = decay_times[decay_times < lmmax]

    ax[0].fill_betweenx(np.linspace(0,1.05,10), lmmin, lmmax, color='darkgray', alpha=0.05)

    frmin = [np.argmin(np.abs(decay_times - lmmax))]
    frmax = [np.argmin(np.abs(decay_times - lmmin))]
    frmin = survive_prob(decay_times, ndcys)[frmin]
    frmax = survive_prob(decay_times, ndcys)[frmax]

    tmin = decay_times[np.argmin(np.abs(survive_prob(decay_times, ndcys) - frmax))]
    tmax = decay_times[np.argmin(np.abs(survive_prob(decay_times, ndcys) - frmin))]
    jfit_times   = lin_fit_times(decay_times, ndcys, tmin, tmax)
    gammas1[tmp] = np.array([temp, -jfit_times[0]])

    col = next(cls)
    ax[0].plot(decay_times, np.exp(get_line(decay_times, *jfit_times)), ls=':', lw=1, color='k')
    ax[0].plot(decay_times, survive_prob(decay_times, ndcys), color=col, ls='-', label=labss)

temprange = np.linspace(tempList[0], tempList[3], 50)
bounce_dat0 = get_best_bounce(tempList[:4], gammas1[:4,1])
ax[1].plot(temprange, exp_bounce(temprange, *bounce_dat0), color='k', label=r'$\rm Best \; Fit$')#label=r'$A \sqrt{\frac{E}{T}} e^{-E/T}$')

temprange = np.linspace(save_teff_from_fld[0], save_teff_from_fld[3], 50)
bounce_dat = get_best_bounce(save_teff_from_fld[:4], gammas1[:4,1])
ax[2].plot(temprange, exp_bounce(temprange, *bounce_dat), color='k', label=r'$\rm Best \; Fit$')#label=r'$A \sqrt{\frac{E}{T_{\rm eff}}} e^{-E/T_{\rm eff}}$')
print(bounce_dat)

cls = cycle(allcolors[:4])
for tmp in range(len(tempList)):
    if tmp>3: break
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])
    col = next(cls)
    ax[1].plot(temp, gammas1[tmp,1], color=col, marker='o', linestyle='none')
    temp2 = save_teff_from_fld[tmp]
    ax[2].plot(temp2, gammas1[tmp,1], color=col, marker='o', linestyle='none')

ax[0].set_yscale('log')
ax[0].set_xlim(0, lmmax)
ax[0].set_ylabel(r'$\rm Pr(survive)$', fontsize=13)
ax[0].set_xlabel(r'$t$')
ax[0].grid(True, ls=':', color='lightgray', alpha=0.5)
ax[0].tick_params(direction='in', which='both', top=True, right=True)

a = tempList[tmpList]
ax[1].set_xticks(a)
a = [round(ii,3) for ii in save_teff_from_fld]
ax[2].set_xticks(a)
for ai, aa in enumerate(ax[1:]):
    aa.set_ylabel(r'$\Gamma / L$', fontsize=13)
    aa.set_xlabel([r'$T$',r'$T_{\rm eff}$'][ai])
    aa.legend(loc=2)
    aa.grid(True, ls=':', color='lightgray', alpha=0.5)
    aa.tick_params(direction='in', which='both', top=True, right=True)
    aa.ticklabel_format(axis='both', style='scientific', scilimits=[-1.,0.])
fig.tight_layout()
plt.show()

/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in sqrt
  """Entry point for launching an IPython kernel.
/cm/shared/apps/python/python37/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: overflow encountered in exp
  """Entry point for launching an IPython kernel.

[0.10948856 0.30947499]
In [34]:
def get_stuffs(nsets):
    partdat = np.zeros((len(tmpList), nsets))
    for ii, average_bubble in enumerate(bubbleList):
        tmp = tmpList[ii]
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        veldata = np.array(np.load(velocities_file(*exp_params)))
        simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]

        pool = set(all_vels)
        means, vars = np.empty(nsets), np.empty(nsets)
        slen = int(round(len(pool) / nsets,1))
        for ss in range(nsets):
            velssec = random.sample(pool, slen)
            pool -= set(velssec)
            velssec = np.array(velssec)
            partdat[tmp,ss] = np.std(velssec)
    return partdat

ecrit, esph, char_vels = np.zeros((3, len(tmpList)))
for ii, average_bubble in enumerate(reversed(bubbleList)):
    tmp = tmpList[ii]
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])

    veldata = np.array(np.load(velocities_file(*exp_params)))
    simvels, all_vels = veldata[:,0], veldata[:,1]
    
    char_vels[tmp] = np.std(all_vels)

    instanton  = instantonList[tmp]
    gradinst   = (np.roll(instanton,-1) - instanton)/dx
    ecrit[tmp] = np.sum(0.5*gradinst**2. + V(instanton, lamb) - V(phieq,lamb))

    esph[tmp]  = encritList[tmp]

diffen = phi0**2.*dx
#sigma_avbub = tempList[tmpList] / (esph * diffen)
#sigma_crit  = tempList[tmpList] / (ecrit * diffen)
#sigma_drate = tempList[tmpList] / bounce_dat[-1]

sigma_avbub = save_teff_from_fld / (esph * diffen)
sigma_crit  = save_teff_from_fld / (ecrit * diffen)
sigma_drate = save_teff_from_fld / bounce_dat[-1]
In [35]:
fig, ax = plt.subplots(1,4, figsize = (8,2), sharey='row')
cols = cycle(allcolors[:4])
save_gauss_fit = []
for ii, average_bubble in enumerate(bubbleList):
    tmp = tmpList[ii]
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])
    labss = r'$T/m = {}$'.format(round(temp/np.sqrt(m2(lamb)),1))#labl(lamb, phi0, temp)

    veldata = np.array(np.load(velocities_file(*exp_params)))
    simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]

    col = next(cols)
    n, bins, patches = ax[ii].hist(all_vels, bins=20, alpha=1, color=col, density=True)#, label=labss)
    n = n[1:]
    bins = bins[1:]
    centers = 0.5*(bins[1:] + bins[:-1])

    xlist = np.linspace(-1.,1.,1000)

    gauss_fit = lambda x, mean, std: np.exp(- (x - mean)**2. / (2.*std**2.)) / np.sqrt(2.*np.pi*std**2.)
    fit_gauss = lambda x, data: sco.curve_fit(gauss_fit, x, data)[0]

    best_gauss = fit_gauss(centers, n)
    ax[ii].plot(xlist, gauss_fit(xlist, *best_gauss), color='k', ls='--')
    save_gauss_fit.append(best_gauss[-1])

    bb = sigma_crit[tmp]**0.5
    ax[ii].plot(xlist, gauss_fit(xlist, 0., bb), color='k', ls='-')

    ax[ii].plot(xlist, gauss_fit(xlist, np.mean(all_vels), np.std(all_vels)), color='k', ls=':')

    a = ax[ii].get_xticks().tolist()[1:-1:]
    a = np.linspace(-1, 1, 5)
    a[::2] = np.array(a[::2], dtype='int')
    ax[ii].set_xticks(a)
    b = [r'${:.1f}$'.format(al) for al in a]
    b[::2] = [r'${:.0f}$'.format(al) for al in a[::2]]
    ax[ii].set_xticklabels(b)
    ax[ii].grid(ls=':', color='darkgray',alpha=0.5)
    ax[ii].set_xlabel(r'$v$')
    ax[ii].set_title(labss, fontsize=10)
    ax[ii].tick_params(direction='in', which='both', top=True, right=True)

a = ax[0].get_yticks().tolist()[:-1:]
a = [round(al,2) for al in a]
ax[0].set_yticks(a)
b = [r'${:.1f}$'.format(al) for al in a]
ax[0].set_yticklabels(b)
ax[0].set_ylabel(r'$\rm PDF$')


plt.savefig('./plots/comparison_vels_PDF.pdf')
plt.show()
In [36]:
fig, ax = plt.subplots(1,1, figsize = (4,2.5))
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), sigma_crit, label=r'$\rm Prediction$')
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), sigma_avbub, label=r'$E_{\rm exp}$')
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), sigma_drate, label=r'$\rm Decay \; Rate$')

for nsets in reversed(np.linspace(15, 20, 1)):
    partdat = get_stuffs(int(nsets))
    ax.errorbar(save_teff_from_fld/np.sqrt(m2(lamb)), np.mean(partdat, axis=-1)**2., yerr=np.std(partdat, axis=-1), \
                alpha=1, lw=1, label=nsets)

ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), char_vels**2., label=r'$\rm Data$')

a = [round(ii,3) for ii in save_teff_from_fld/np.sqrt(m2(lamb))]
ax.set_xticks(a)
ax.grid(True, ls=':', color='lightgray', alpha=0.5)
ax.tick_params(direction='in', which='both', top=True, right=True)
ax.ticklabel_format(axis='both', style='scientific', scilimits=[-1.,0.])
ax.legend(bbox_to_anchor=(1, 1))

ax.set_xlabel(r'$T_{\rm eff}/m$')
ax.set_ylabel(r'$\left<v^2\right>$')
plt.tight_layout()
plt.show()
In [37]:
fig, ax = plt.subplots(1,1, figsize = (5.5,3))

ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), sigma_crit, label=r'$\rm Prediction$')

partdat = get_stuffs(15)
ax.errorbar(save_teff_from_fld/np.sqrt(m2(lamb)), char_vels**2., \
            yerr=np.std(partdat, axis=-1), alpha=1, lw=1, label=r'$\rm Variance \; of \; Data$')
ax.errorbar(save_teff_from_fld/np.sqrt(m2(lamb)), np.array(save_gauss_fit)**2., \
            yerr=np.std(partdat, axis=-1), alpha=1, lw=1, label=r'$\rm Gaussian \; Best \; Fit$')

a = [round(ii,3) for ii in save_teff_from_fld/np.sqrt(m2(lamb))]
ax.set_xticks(a)
ax.grid(True, ls=':', color='lightgray', alpha=0.5)
ax.tick_params(direction='in', which='both', top=True, right=True)
ax.ticklabel_format(axis='x', style='scientific', scilimits=[-1.,0.])
ax.ticklabel_format(axis='y', style='scientific', scilimits=[-10.,10.])
a = ax.get_yticks().tolist()[1::2]
a = [round(al,2) for al in a]
ax.set_yticks(a)
a = [r'${:.2f}$'.format(al) for al in a]
ax.set_yticklabels(a)

leg = ax.legend(bbox_to_anchor=(1.5, 0.6))._legend_box.align='left'

ax.set_xlabel(r'$T_{\rm eff}/m$')
ax.set_ylabel(r'$\left<v^2\right>$')
plt.tight_layout()
plt.savefig('./plots/comparison_vcom_pred.pdf')
plt.show()
In [38]:
#velsq_exp = lambda edt: 1. - edt - edt**2.* np.exp(edt) * scp.special.expi(-edt)
#velsq_exp_explicit = lambda en, tm: 1. - en/tm - (en/tm)**2.* np.exp(en/tm) * scp.special.expi(-en/tm)

#def f_inv_velsq_exp_explicit(x, ii):
#    return abs(velsq_exp_explicit(x, tempList[tmpList][ii]) - char_vels[ii])

#enfromvels = np.zeros((len(char_vels)))
#for ii in range(len(char_vels)):
#    out = sco.minimize_scalar(f_inv_velsq_exp_explicit, bounds=((0,None)), args=(ii))
#    enfromvels[ii] = out.x
In [39]:
# Plot to compare all energies 

fig, ax = plt.subplots(1, 1, figsize = (5.5,3.))
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), 1./sigma_crit, label=r'$\rm Prediction$')
ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), 1./sigma_avbub, label=r'$\rm Critical \; Bubble$')

ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), 1./sigma_drate, label=r'$\rm Decay \; Rate$')

ax.plot(save_teff_from_fld/np.sqrt(m2(lamb)), 1./np.array(save_gauss_fit)**2., label=r'$\rm Velocity \; Distribution$')

a = [round(ii,3) for ii in save_teff_from_fld/np.sqrt(m2(lamb))]
ax.set_xticks(a)
ax.ticklabel_format(axis='x', style='scientific', scilimits=[-1.,0.])
ax.ticklabel_format(axis='y', style='scientific', scilimits=[-10.,10.])

ax.set_ylabel(r'$E/T_{\rm eff}$')
ax.set_xlabel(r'$T_{\rm eff}/m$')

ax.legend(ncol=1, loc='best', frameon=False, bbox_to_anchor=(1, 1))._legend_box.align='left'
ax.grid(which='both', color='darkgray', ls=':', alpha=0.4)
ax.tick_params(direction='in', which='both', top=True, right=True)

plt.tight_layout()
plt.savefig('./plots/critenergy_comparison.pdf')
plt.show()
In [40]:
for ii, average_bubble in enumerate(bubbleList):
    tmp = tmpList[ii]
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])

    decay_times = np.load(decay_times_file(*exp_params, minSim, maxSim, nTimeMAX))
    print('N1 = ', len(decay_times))

    minDecTime  = nLat*2//3
    alltimes    = decay_times[:,1]
    simList2Do  = decay_times[alltimes>=minDecTime, 0]
    n2Do        = len(simList2Do)
    print('N2 = ', n2Do)

    veldata = np.array(np.load(velocities_file(*exp_params)))
    simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]
    
    print('N3 = ', len(all_vels))
 #   print(temp, len(all_vels), np.count_nonzero(all_vels > 0.)/len(all_vels))
  #  print(temp, np.count_nonzero(np.abs(all_vels) > 0.5)/len(all_vels))

N1 =  1003
N2 =  665
N3 =  654
N1 =  2192
N2 =  1199
N3 =  1191
N1 =  3165
N2 =  1280
N3 =  1266
N1 =  3711
N2 =  882
N3 =  854
In [41]:
if False:
    for tmp in range(len(tempList)):
        if tmp!=0: continue

        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        # 0 for field; 1 for momentum 
        find = 0
        # modes to plot:
        aa, bb = 1, knyq-1
        labs = labl(lamb, phi0, temp)
        psth = pspec(lamb, phi0, temp)[aa:bb]

        ALL_powspec1 = np.load(powspec_tlist_file(*exp_params, minSim, 2000))
        ALL_powspec2 = np.load(powspec_tlist_file(*exp_params, 2000, maxSim))
        tlist, PSfld1, PSfld2 = ALL_powspec1[0], ALL_powspec1[1], ALL_powspec2[1]

        fig, ax = plt.subplots(1,1, figsize = (8,4))
        camera = Camera(fig)
        for tind, tt in enumerate(tlist):
            avPSfld1 = np.nanmean(PSfld1[:, find, tind, :], axis=0)
            avPSfld2 = np.nanmean(PSfld2[:, find, tind, :], axis=0)
            avPSfld  = np.nanmean(np.array([avPSfld1[:], avPSfld2[:]]), axis=0)
            slice = plt.plot(klist[aa:bb], avPSfld[aa:bb]/psth, ls='-', color='purple')

            ax.set_xscale('log')
            ax.set_yscale('log')
            ax.set_ylabel(r'$P(k)$')
            ax.set_xlabel(r'$k$')
            ax.legend(slice, [f't = {tt}'], loc=1, title=labs)
            camera.snap()
        animation = camera.animate(interval = 0.05);
        animation.save('./plots/animation_cut_PS'+batch_params(*exp_params)+'.gif', writer = 'imagemagick')
In [42]:
fig, ax = plt.subplots(1,1, figsize = (5,2.5))
for tmp in range(len(tempList)):
    phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
    exp_params = np.asarray([nLat, lamb, phi0, temp])

    aa, bb = 1, knyq-1
    plt.plot(klist[aa:bb], pspec(lamb, phi0, temp)[aa:bb], ls=lsl(tmp), label = labl(lamb, phi0, temp)) # th pow spec
plt.xscale('log')
plt.legend(); plt.show()
In [43]:
if True:
    choose = random.sample(np.arange(4000).tolist(), 4)
    for sim in choose:
        fig, ax = plt.subplots(1,1, figsize = (4,2.5))
        for tmp in range(len(tempList)):
            phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
            exp_params = np.asarray([nLat, lamb, phi0, temp])

            if sim < 2000: minSim, maxSim = 0, 2000
            else: minSim, maxSim, sim = 2000, 4000, sim-2000

            ALL_toten = np.load(toten_tlist_file(*exp_params, minSim, maxSim))
            tlist, energy = ALL_toten[0], ALL_toten[1][sim]

            try: nnrg = np.argwhere(np.isnan(energy)).flatten()[0]
            except: nnrg = len(tlist)+1

            tcut, tencut = tlist[:nnrg], energy[:nnrg]
            tencut = dx * phi0**2. * (tencut - tencut[0]) / tencut[0]
            plt.plot(tcut, tencut, label=labl(lamb, phi0, temp))
        plt.legend(title='Sim='+str(sim), bbox_to_anchor=(1.2,1.))
        plt.show()
In [44]:
if True:
    for tmp in range(len(tempList)):
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])
        labs = labl(lamb, phi0, temp)
        
        minSim = 2000

        ALL_toten = np.load(toten_tlist_file(*exp_params, minSim, maxSim))
        ALL_emt   = np.load(emt_tlist_file(*exp_params, minSim, maxSim))

        tlist, enfld, emtfld  = ALL_toten[0], ALL_toten[1], ALL_emt[1]

        fig, ax = plt.subplots(1,1, figsize = (4,2.5))
        for sim, (momentum, energy) in enumerate(zip(emtfld, enfld)):
            if sim!=0: continue

            try: nnrg = np.argwhere(np.isnan(energy)).flatten()[0]
            except: nnrg = len(tlist)+1

            tcut, emtcut, tencut = tlist[:nnrg], momentum[:nnrg], energy[:nnrg]
            emtcut = (emtcut - emtcut[0]) / emtcut[0]
            tencut = (tencut - tencut[0]) / tencut[0]

            ax.plot(tcut, tencut, label=(r'$T^{00}$' if sim==0 else None))
            ax.plot(tcut, emtcut, label=(r'$T^{0x}$' if sim==0 else None))

        ax.set_xlabel(r'$\phi_0^{-1} \sqrt{V_0} \; t$')
        ax.set_ylabel(r'$T^{00}, \; T^{0x}$')
        ax.legend(title=labs, bbox_to_anchor=(1.,1.))
        ax.grid(ls=':', color='darkgray', alpha=0.3)
        ax.set_title(labs)
        fig.show()
In [45]:
if True:
    fig, ax = plt.subplots(1, 1, figsize=(3.5, 2.5))
    allcs = cycle(allcolors)
    for tmp in reversed(range(len(tempList))):
        if tmp>3: continue

        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])
        labss = r'${}$'.format(round(temp/np.sqrt(m2(lamb)),1))#labl(lamb, phi0, temp)
        col = next(allcs)

        veldata = np.array(np.load(velocities_file(*exp_params)))
        simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]

        ALL_emt1 = np.load(emt_tlist_file(*exp_params, minSim, 2000))
        ALL_emt2 = np.load(emt_tlist_file(*exp_params, 2000, maxSim))
        emtfld = np.concatenate((ALL_emt1[1], ALL_emt2[1]), axis=0)
        initemt    = emtfld[simvels,0]

        #emtfld = np.load(stdemt0_tlist_file(*exp_params, minSim, maxSim))
        #initemt    = emtfld[simvels]

        ax.plot(all_vels, -initemt, marker='o', ms=2, linestyle='None', label=labss, color=col)
 
        f = lambda x, A: x/A
        pbf, pcov = sco.curve_fit(f, all_vels, -initemt) # your data x, y to fit
        ax.plot(all_vels, f(all_vels, *pbf), color='k', linewidth=1)
    
    ax.set_xlim(-1,1)
  #  ax.set_ylim(-1,1)
    ax.set_xlabel(r'$v_{\rm COM}$')
#    ax.set_ylabel(r'$\int \mathrm{d} r \, \Pi \, \partial_r \varphi$')
    ax.set_ylabel(r'$\bar{P}$')
    ax.grid(ls=':', color='darkgray', alpha=0.3)

    leg = ax.legend(title=r'$T/m:$', ncol=1, loc='best', fontsize=10, bbox_to_anchor=(1, 0.8))

    if True:
        a = ax.get_yticks().tolist()[2:-1:2]
        a = [round(al,2) for al in a]
        ax.set_yticks(a)
     #   a[::2] = [int(al) for al in a[::2]]
     #   b = [r'${:.1f}$'.format(al) for al in a]
     #   b[::2] = [r'${:.0f}$'.format(al) for al in a[::2]]
     #   ax.set_yticklabels(b)

    a = ax.get_xticks().tolist()
    a = [round(al,2) for al in a]
    ax.set_xticks(a)
    a[::2] = [int(al) for al in a[::2]]
    b = [r'${:.1f}$'.format(al) for al in a]
    b[::2] = [r'${:.0f}$'.format(al) for al in a[::2]]
    ax.set_xticklabels(b)

    plt.savefig('./plots/initial_emt_vs_vel.pdf')
    fig.show()
In [46]:
if True:
    fig, ax = plt.subplots(1, 1, figsize=(5, 4))
    allcs = cycle(allcolors)
    for tmp in reversed(range(len(tempList))):
        if tmp>3: continue

        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])
        labss = r'${}$'.format(round(temp/np.sqrt(m2(lamb)),1))#labl(lamb, phi0, temp)
        col = next(allcs)

        veldata = np.array(np.load(velocities_file(*exp_params)))
        simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]

        ALL_emt1 = np.load(emt_tlist_file(*exp_params, minSim, 2000))
        ALL_emt2 = np.load(emt_tlist_file(*exp_params, 2000, maxSim))
        emtfld = np.concatenate((ALL_emt1[1], ALL_emt2[1]), axis=0)
        initemt    = emtfld[simvels,0]

        # fit a straight line
        jfit_times = np.polyfit(all_vels, -initemt, deg=1)
        print(jfit_times)
        vlist = np.linspace(0., 1, 100)
        ax.plot(vlist, get_line(vlist, *jfit_times), '-', label=labss, color=col)

    ax.legend(title=r'$T/m {\rm \; best \; fit}$', loc='best', ncol=2)
    ax.set_xlabel(r'$v_{\rm COM}$')
    ax.set_ylabel(r'$\bar{P}$')
    ax.grid(ls=':', color='darkgray', alpha=0.3)
    fig.show()

[ 0.18583269 -0.00668639]
[ 0.18248509 -0.00590275]
[0.18450292 0.00386879]
[0.18669835 0.00104476]
In [47]:
if False:
   # for tmp in reversed(range(len(tempList))):
    for tmp in reversed(range(len(tempList))):
        if tmp>3: continue
        
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])

        path = decay_times_file(*exp_params, minSim, maxSim, nTimeMAX)
        if os.path.exists(path):
            print(path)
            decay_times = np.load(path)

            minDecTime = nLat*2//3
            alltimes   = decay_times[:,1]
            simList2Do = decay_times[alltimes>=minDecTime, 0]

            all_vels = []
            for sim in simList2Do:
                path2RESTsim = rest_sim_location(*exp_params, sim)
                if os.path.exists(path2RESTsim):

                    sim, bubble, totbeta = np.load(path2RESTsim)
                    all_vels.append(np.array([sim, totbeta]))

        np.save(velocities_file(*exp_params), all_vels)
        print(len(all_vels), len(simList2Do))
        print('Done!')
In [48]:
if True:
    fig2, ax = plt.subplots(1, 1, figsize=(6, 5))
    plt.grid(True, ls=':', alpha=0.5)
    for tmp in reversed(range(len(tempList))):
        if tmp>3: continue
    
        phi0, temp, lamb, sigmafld, minSim, maxSim, right_Vmax = get_model(tmp)
        exp_params = np.asarray([nLat, lamb, phi0, temp])
        labss = r'${}$'.format(round(temp/np.sqrt(m2(lamb)),1))#labl(lamb, phi0, temp)

        veldata = np.load(velocities_file(*exp_params))
     
        simvels, all_vels = veldata[:,0].astype(int), veldata[:,1]
        print(temp, len(simvels))
        plt.hist(all_vels, bins=12, label=labss, density=True)
    plt.xlim((-1,1))
    plt.legend(title=r'$T/m$', ncol=1, loc='best')
    plt.xlabel(r'$v$')
    plt.legend()
    plt.show()

0.12 854
0.11 1266
0.1 1191
0.09 654
In [ ]:
In [ ]: